function varargout = PsychGenPlot(S,varargin)
% PsychGenPlot(S);
% PsychGenPlot(S,'propertyname',propertyvalue);
% HR = PsychGenPlot(S,...);
% [HR,FAR] = PsychGenPlot(S,...);
% [HR,FAR,cHR] = PsychGenPlot(S,...);
% [HR,FAR,cHR,Sens] = PsychGenPlot(S,...);
% 
% S is the schedule structure built using the ScheduleMaker GUI.
%
% Possible property input pairs:
%   'Parent',    target axes (default = gca)
%   'dependant', char string corresponding to a value from S.readparams 
%                   ex: 'sound_level'
%   'condition', char string corresponding to a value from S.readparams 
%                   ex: 'carrier_frequency'
%   'xscale',    either 'linear' or 'log' (default)
% 
% S can also contain a plotcfg sub-structure to set the ploting parameters
% for the schedule S.  Any values explicitly defined during a call to this
% function using property pairs (see abov) will have priority over whatever
% was predifined in the plotcfg sub-structure.  Note that only the 'Parent'
% property cannot be set in the plotcfg sub-structure.  Field names in the
% plotcfg sub-structure should be the same as the property names listed
% above.
% 
% 
% Outputs:
%   HR:     hit rate
%   FAR:    False-alarm rate
%   cHR:    False-alarm corrected hit rate (Heffner & Heffner, 1988)
%   Sens:   50% sensitivity based on cHR
% 
% DJS (c) 2011

ax  = gca;
xsc = 'log';
showlegend = true;

% check for and process plotcfg substructure in S
if isfield(S,'plotcfg')
    fn = fieldnames(S.plotcfg);
    for i = 1:length(fn)
        switch lower(fn{i})
            case 'dependant'
                dvar = S.plotcfg.(fn{i});
            case 'condition'
                cvar = S.plotcfg.(fn{i});
            case 'showlegend'
                showlegend = S.plotcfg.(fn{i});
            case 'xscale'
                xsc  = S.plotcfg.(fn{i});
        end
    end
end

% check property pairs (trumps S.plotcfg properties)
for i = 1:2:nargin-1
    switch lower(varargin{i})
        case 'parent'
            ax   = varargin{i+1};
        case 'dependant'
            dvar = varargin{i+1};
        case 'condition'
            cvar = varargin{i+1};
        case 'showlegend'
            showlegend = varargin{i+1};
        case 'xscale'
            xsc  = varargin{i+1};
    end
end



if isempty(ax), ax = gca; end

h = 2;      % hits
% m = 1;      % misses
f = 3;      % false alarm

rv = S.response_vals;
if ~any(S.response_vals(end,2:end)) || all(isnan(S.response_vals(end,2:end)))
    rv(end,:) = []; % last response value is invalid
end

% find trial types
t = strfind(S.readparams,'trial_type');
xtt = findincell(t,1) + 1;
catch_trials = rv(:,xtt) == 1;

% find response codes
t = strfind(S.readparams,'response_code');
xrc = findincell(t,1) + 1;

% find condition
t = strfind(S.readparams,cvar);
xcvar = findincell(t,1) + 1;
ucvar = unique(rv(~catch_trials,xcvar));

% find dependant variable
t = strfind(S.readparams,dvar);
xdvar = findincell(t,1) + 1;
udvar = unique(rv(~catch_trials,xdvar));

% compute hit rate for each Frequency & DB level
hr = zeros(length(udvar),length(ucvar));
for i = 1:length(ucvar)
    for j = 1:length(udvar)
        ind   = rv(:,xcvar) == ucvar(i) & rv(:,xdvar) == udvar(j);
        hind  = ismember(rv(ind,xrc),h);
        hr(j,i) = sum(hind)/sum(ind);
    end
end
hr(isnan(hr)) = 0;

% compute false-alarm rate
fa = sum(rv(catch_trials,xrc) == f) / sum(catch_trials);

% plot data
cla(ax);

% plot(ax,udvar,hr,'LineWidth',2);
plot(ax,udvar,hr,'o-','MarkerSize',3);

hold(ax,'on');

% compute corrected hit rate and estimate sensitivity
Sens = nan(1,size(hr,2)); cHR = nan(sum(udvar>0),size(hr,2));
if size(hr,1) > 2 && any(any(hr))
    for i = 1:size(hr,2)
        cHR(:,i) = hr(udvar>0,i) - hr(udvar>0,i) * fa; % corrected hit rate (Heffner & Heffner, 1988)
        Sens(i)  = interp1q(cHR(:,i),udvar(udvar>0),0.5);  % estimate sensitivity
    end
end

ylim(ax,[0 1]);

set(ax,'xscale',xsc);

plot(ax,xlim(ax),[fa fa],':r');

hold(ax,'off');

xlabel(ax,dvar);
ylabel(ax,'Hit Rate');

if showlegend
    h = legend(ax,num2str(ucvar/1000,'%1.0fkHz'),'Location','SouthEast');
    set(h,'Box','off');
    set(h,'Color','none');
end

% outputs
varargout{1} = hr;
varargout{2} = fa;
varargout{3} = cHR;
varargout{4} = Sens;




